Symmetry breaking and criticality in tensor-product states 
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We discuss variationally optimized matrix-product states for the transverse-field Ising chain, using 
D x D matrices with small D £ {2 — 10}. For finite system size N there are energy minimums for 
symmetric as well as symmetry-broken states, which cross each other at a field value h c (N,D); 
thus the transition is first-order. A continuous transition develops as N — > oo. The asymptotic 
critical behavior is then always of mean- field type (the magnetization exponent (3 = 1/2), but a 
window of field strengths where true Ising scaling holds (/? = 1/8) emerges with increasing D. We 
also demonstrate asymptotic mean-field behavior for infinite-size two-dimensional tensor-product 
(iPEPS) states with small tensors. 

PACS numbers: 75.10.Jm, 75.40.Mg, 75.40.Cx, 05.30.Rt 



Methods based on matrix-product states (MPSs) [HQ 
have become the primary computational tools for studies 
of static as well as dynamic properties of one-dimensional 
quantum many-body systems [3(. Key steps in the de- 
velopment of these methods were White's density matrix 
renormalization group (DMRG) 0,1, the demonstra- 
tion by Ostlund and Romer of its connection to MPSs 
0, and later important insights from the held of quan- 
tum information theory. In particular, the concept of 
entanglement entropy (the area law) both explains the 
success of the approach in one dimension and its failure 
(violation of the area law) in higher dimensions j3~§- 
The formulation of computational methods directly in 
terms of MPSs also led to a framework for efficient op- 
timization of these states independently of the DMRG 
method 0-U|, and to a long-sought way of computing 
time evolution The MPS approach also has a nat- 
ural extension to higher dimensions which does obey the 
area law 0] — tensor-product states; also referred to as 
projected-entangled-pair-states (PEPSs) 

In spite of numerous successful applications of MPS- 
based methods, some fundamental aspects of this class of 
quantum states have not yet been studied in detail. It is 
well known that the finite size D of the D x D matrices 
(the elements of which are the variational parameters) 
imposes a finite correlation length, and recently it has 
been recognized that scaling in D for infinite system size 
N can be carried out as an alternative to finite-size scal- 
ing [l5| (i.e., D and N can be considered as different 
but equally valid ways to regularize the calculations) . As 
in mean- field theory (which corresponds to D — 1), an 
MPS can break symmetries of the hamiltonian at a phase 
transition. Exactly how the critical behavior of the or- 
der parameter (the true scaling exponent f3) emerges as 
a function of N and D has not been studied systemat- 
ically, however. This may be partially due to technical 
challenges in properly optimizing an MPS close to a phase 
transition. Such issues are present also for the PEPS ap- 



proach in two dimensions. Order-parameter curves often 
exhibit rounding [171 ] , that may appear due to incomplete 
convergence, approximations made |15|, or due to exter- 
nal fields included to stabilize the calculation [l8[. Nev- 
ertheless, the behavior slightly away from the transition 
can be well described by the expected critical exponent 
The question remains whether this is the 
true critical behavior of the MPS or PEPS variational 
ansatz with finite D, or whether there could eventually 
be a cross-over to a different asymptotic form. 

In this Letter we study the asymptotic critical behav- 
ior by using numerically stable high-precision optimiza- 
tion methods for small D, for both finite and infinite N. 
Using the transverse-field Ising model as a demonstra- 
tion, we show that access to the true critical behavior of 
an MPS requires very high numerical precision; in some 
cases higher than the double-precision (64-bit) floating 
point arithmetic normally used. Optimizing the states 
to the required precision, we show that the asymptotic 
critical behavior of the order parameter is always mean- 
field like ((3 — 1/2). The true universal exponent for the 
model (/? = 1/8) emerges in a window which approaches 
the critical point as D increases. We also shows results 
in two dimensions for an infinite-size PEPS (iPEPS), 
optimized using a recently developed numerically stable 
scheme [23]. Also here we find /3 = 1/2 asymptotically. 

First, consider the simplest kind of MPS for a periodic, 
translationally invariant S = 1/2 spin chain; 
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where a? = ±1 and A(±l) are two hermitian D x D ma- 
trices. As illustrated in Fig. [TJ the normalization of this 
state can be expressed as the contraction of a network 
of 3-index tensors A a b(cr), where a = ±1 is the physi- 
cal index. By contracting over the physical indices first, 
matrices B of size D 2 x D 2 are obtained; 

Ba = A ab (+l)A* cd (+l) + A ab {-l)A: d (-l), (2) 
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FIG. 1: (Color online) (a) The norm of an MPS expressed as 
the contraction of a tensor network. Carrying out the summa- 
tions over the spin indices first (vertical bonds), as indicated 
in (b), gives a simple trace of a product of matrices of size 
D 2 x D 2 (with a possible labeling of the elements indicated). 



where i — a + (c — 1)D and j = i 
normalization is then simply 

= Tr{B N } 



(d - 1)D. The 



(3) 



Expectation values can be computed in a very similar 
way, with some of the B matrices in the product replaced 
by the matrix obtained as in @ but with the operator 



in question first acting on the physical index 
instance, the magnetization m is given by 

Tt{MB n - 1 } 
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where the matrix M is 



M« = A ab (+l)A* cd (+l) - A ab {-l)A* cd {-l). (5) 

The generalization to expectation values of products of 
two or more operators is straight-forward. 

The matrix B is exactly analogous to the transfer ma- 
trix in classical statistical mechanics. With U the uni- 
tary matrix that diagonalizes B, giving its eigenvalues 
Ai, . . . , X]j2, the magnetization can be written as 



^[U^MUjuX 
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As in the transfer-matrix approach, the N — > oo limit can 
be taken by keeping only the leading eigenvalue; assumed 
here to be Ai . The magnetization is then 
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(7) 



where v\ is the eigenvector of B corresponding to Ai. 

Given a hamiltonian H, the problem is how to find the 
matrices A(±l), of given size D, that best reproduce the 
ground state. This can be formulated as a variational 
problem; to minimize the energy E = (^f\H\^f). Several 
different optimization methods have been developed. For 
finite N, the translational invariance is typically broken 
as a series of local optimizations are carried out, sweep- 
ing back and forth through an open chain Q (similar to 
DMRG calculations [1, Q). In a periodic chain, where 
the calculation is more demanding, uniformity is gradu- 
ally restored as the matrices converge. For N = oo, the 
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FIG. 2: (Color online) Stochastic energy minimization with 
10 4 updates per step for a D = 4 MPS at h/J = 1.01432, 
using 64-bit floating point arithmetic. The relative energy 
and magnetization errors are defined as Ae = (E — Eoa)/Eoo, 
A m = \\m\ ~ ImooH/lmool, where the subscript oo refers to 
results converged at the 128-bit level. 

most efficient approach is Vidal's time evolving block dec- 
imation (TEBD) scheme (T^j, where the ground state is 
projected out in the limit of long imaginary time [3, E| > 
starting from an initial (e.g., random) state. Similar 
methods have also been developed for two-dimensional 
iPEPSs, where expectation values cannot simply be ex- 
pressed in eigenvalue forms such as (JS]), but where good 
approximations to the contractions can still be defined 
and evaluated using TEBD-like methods [l6j |. 

Here we investigate symmetry breaking and critical 
scaling of the order parameter in the transverse- field Ising 
model. In one dimension the hamiltonian is 



H = 



N 
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with periodic boundary conditions. This model is exactly 
solvable 2lj and has a paramagnetic-magnetic (m ^ 0) 
transition at h c /J =1. In two dimension, the critical 
point has been determined using quantum Monte Carlo 
calculations, giving h c j J w 3.044 22]. Single-spin mean- 
field theory (D — 1, N — oo) gives h c /J = 2 and 4 in 
one and two dimensions, respectively, and the mean- field 
form of the magnetization is m ~ (h c — h)@ for h < h c , 
with j3 = 1/2. The exact critical exponent is /3 = 1/8 in 
one dimension and /3 rs 0.325 in two dimensions. 

Considering first MPSs, we optimize the A matrices 
(which we take as real and symmetric) using two dif- 
ferent stochastic schemes; one using derivatives and one 
using only the energy. While the convergence is very slow 
for large D compared to state-of-the-art TEBD ll|, the 
methods do not rely on any approximations and are nu- 
merically stable. With stochastic updates, we can avoid 
potential local minimums in a complex energy landscape. 
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FIG. 3: (Color online) Magnetization curves for D = 2 and 
different system sizes. The solid curve is for N = oo. The 
inset shows two almost degenerate energy minimums for TV = 
12, which cross each other at the transition. 

The derivative-based method is the one used in [24]], but 
with exactly computed energies and derivatives for finite 
N. For N = oo, we instead use a brute-force scheme 
with completely random simultaneous updates of all the 
matrix elements (but keeping the matrices symmetric); 
^ab(c) — > A a b(a) + <5[l/2 — r a b(c)], with uniformly dis- 
tributed random numbers r a b(a) G [0,1). An update is 
accepted only if the energy decreases, and then the matri- 
ces are normalized so that the largest |A a b(er)| = 1. One 
step of this procedure typically involves n ~ 10 3 — 10 4 
trials. If the acceptance rate is below 10% we reduce S by 
dividing by, e.g., 1.1. To ensure full convergence, when 8 
has reached the limit where the updates no longer can in- 
fluence the energy (within the numerical precision), it is 
reset to a larger value and the process is repeated (several 
times, until no updates are accepted). 

Fig. [5] illustrates the brute- force procedure for a D = 4 
MPS optimized at h/J = 1.01432. The evolution of 
the errors of the energy and the magnetization is shown, 
along with S and the acceptance rate. In this case the ac- 
ceptance rate was always below 10%, and 6 therefore de- 
creases after each step. This calculation was carried out 
using standard 64-bit floating-point arithmetic, which is 
reflected in the convergence of the energy to within a rel- 
ative error of s=s 10 -15 . The computation was continued 
with 128-bit arithmetic, until the energy was converged 
to w 10~ 25 . The errors graphed in the figure are with re- 
spect to this second optimization. The 64-bit optimiza- 
tion took only a few minutes, whereas the subsequent 
128-bit run took many hours. The computational effort 
increases very rapidly with D, and we have only carried 
out systematic studies up to D = 10 (for which some 
points required several weeks of CPU time) [23| . 

Note that while the energy in Fig. [2] has converged to 
full 64-bit precision, the relative magnetization error is 
much larger, A m w 10~ 3 . Using 128-bit arithmetic gives 
m = 0.031814167 (where all digits shown are converged). 
It is well known that the energy in MPS and DMRG cal- 
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FIG. 4: (Color online) Demonstration of asymptotic MPS 
mean-field behavior and scaling cross-over in one dimension. 
The D-dependent critical fields are: h c /J = 1.0717967 (D — 
2), 1.0143343 (D = 4), 1.0063523 (D = 6) 1.0021654 (D = 
10). The lines have slopes (3 = 1/8 and 1/2. 

culations converges much faster than other quantities [j| , 
but N — oo results close to the critical point appear to be 
even more sensitive to extremely small energy variations 
than had been previously anticipated. When trying to 
extract the asymptotic critical behavior of m, the prob- 
lem is accentuated by the fact that it is the relative, not 
absolute error that is relevant. All results to be discussed 
below have been converged to the level required for a re- 
liable scaling analysis. 

As shown in Fig. [31 for finite N the phase transition oc- 
curs with sharp magnetization jumps for small N, which 
become less pronounced as N increases and the transi- 
tion moves toward higher h. The curves converge toward 
the continuous transition obtained in the infinite- /V cal- 
culation. The first-order behavior can be traced to the 
presence of two energy minimums (shown in Fig. [3J for 
N = 12), which we can track using steepest-decent op- 
timizations starting from large and small h (changing h 
slowly) . The diminishing discontinuity with increasing N 
implies that the minimums move closer to each other in 
parameter space, coinciding at h c for N = oo. For fixed 
finite JV, the discontinuous jumps move toward h = 
with increasing D, reflecting the fact that when D — » oo 
an MPS can reproduce the exact spin-inversion symmet- 
ric (m = 0) ground state of a finite chain. 

For N = oo and any D, the optimal state is symmetry- 
broken below some h c (D), with h c (D)/J 1 as D — > oo. 
The D dependence is not smooth, as has been pointed 
out before [THJ . Here we focus on the behavior of m 
for h — > h c (D). Thanks to our high-precision data, we 
can extract h c {D) reliably using a power-law assumption; 
m oc (h c — h)P for < m <C 1. This always gives (3 w 
0.50 for the best fit, suggesting that the MPS procedure 
leads to mean-field behavior for any finite D. As shown 
in Fig. 01 the true critical behavior (/? = 1/8) emerges 
within a window of h- values with increasing D, with the 
cross-over to j3 — 1/2 gradually moving toward h c . 
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FIG. 5: (Color online) Field dependence of the magnetiza- 
tion computed with a D = 2 iPEPS in two dimensions. The 
critical field is h c /J — 3.1041. 



It is perhaps not surprising, after all, that a finite- 
ly MPS cannot reproduce a non-trivial critical expo- 
nent asymptotically, because the correlation length is fi- 
nite. Criticality (which can be non-mean-field) in a one- 
dimensional classical Ising model requires long-range in- 
teractions [2HI I and the partition function then does not 
correspond to an MPS with finite D. It has also been 
proved that a finite-Z? MPS can be renormalized to a 
product state (2(|. It is, however, remarkable that the 
system is so sensitive to incomplete optimization that 
the asymptotic mean-field behavior of the order param- 
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eter had not been noted in previous studies 

We now turn to the two dimensional iPEPS. Non- 
trivial criticality has been anticipated in this case, even 
for finite D 0], because partition functions of classical 
models with critical points can be written as tensor prod- 
ucts [13J. Magnetization curves closely following the ex- 
pected pow er-law with j3 ~ 0.325 have been reported 
[HI [l7l Il9|. but the calculations are not very accurate 
close to the critical point. We have used an improved 
iPEPS contraction scheme which is less affected by 
approximations. Fig. [5] shows transverse-field Ising re- 
sults for D = 2. An asymptotic mean-field behavior is 
seen unambiguously, and further away from the critical 
point there is again a cross-over to a behavior matching 
closer the true /?. However, for D = 2 the cross-over takes 
place where m is already large, w 0.5, and this is not ac- 
tual critical behavior. Unless D is much larger, there will 
be no clear-cut scaling with the correct exponent. 

Our study shows that great care has to be taken when 
extracting critical scaling forms from the order parame- 
ter in MPS and PEPS calculations. Asymptotic mean- 
field behavior should be expected, not just for the Ising 
models considered here, but at symmetry-breaking tran- 
sitions in general. This information helps to accurately 
locate the critical point for small D. To extract the true 
critical behavior, it is necessary to carefully examine the 
behavior for increasing D. 
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